Simulation of phase transitions in highly asymmetric fluid mixtures 
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We present a novel method for the accurate numerical determination of the phase behavior of fluid 
mixtures having large particle size asymmetries. By incorporating the recently developed geometric 
cluster algorithm within a restricted Gibbs ensemble, we are able to probe directly the density and 
concentration fluctuations that drive phase transitions, but that are inaccessible to conventional 
simulation algorithms. We develop a finite-size scaling theory that relates these density fluctuations 
to those of the grand- canonical ensemble, thereby enabling accurate location of critical points and 
coexistence curves of multicomponent fluids. Several illustrative examples are presented. 
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The vast majority of commercially relevant fluids are 
multicomponent mixtures. An understanding of the 
phase behavior of these systems is of paramount impor- 
tance for applications, and also a matter of great funda- 
mental interest . With the advent of powerful comput- 
ers, various computational techniques have been devised 
to directly determine fluid phase behavior 2]. However, 
these methods are all restricted to fluids in which the 
various components have similar sizes, whereas impor- 
tant phenomena occur in highly size-asymmetric multi- 
component fluids such as colloidal dispersions, colloid- 
nanoparticle mixtures, and polymer solutions 

The computational bottleneck for existing simulation 
methods arises from the fact that they cannot simulta- 
neously relax a fluid system on disparate length scales. 
Specifically, for large size ratios, the big particles be- 
come 'jammed' by the smaller ones. Recently, however, 
building on earlier work 0|, a geometric cluster Monte 
Carlo algorithm (GCA) was proposed that facili- 
tates rejection- free simulations of highly size-asymmetric 
mixtures via large-scale collective updates which move 
whole groups of particles in a single step. Although this 
method has been successfully applied to problems relat- 
ing to colloidal stabilization Ig j in which size asym- 
metries span several orders of magnitude, it is incapable 
of dealing with the density and concentration fluctua- 
tions associated with phase separation, since it inherently 
operates in a canonical ensemble in which the particle 
number and volume are fixed. This limitation cannot 
be overcome by incorporating cluster moves into a stan- 
dard grand-canonical (GC) or constant-iVpT ensemble, 
because this does not address the underlying shortcom- 
ings of these ensembles with respect to sampling of den- 
sity fiuctuations in asymmetric mixtures. 

It is the purpose of this Letter to introduce a method 
that overcomes these problems. This is achieved by em- 
beddiiig cluster moves in a variant of the Gibbs ensem- 
ble [Sl llfll , in such as way that they couple to the density 
fluctuations, resulting in efficient exploration of config- 
uration space. To exploit this approach we present a 



finite-size scaling theory that permits the determination 
of the critical point and the phase boundary. As an il- 
lustration, we apply the method to study liquid-vapor 
coexistence in asymmetric binary mixtures, for which we 
show that the presence of even small quantities of small- 
particle additives can strongly affect the location of the 
critical point. Furthermore, depending on the nature of 
the interaction of the additive with the fiuid particles, the 
critical temperature can either be enhanced or depressed. 

To enable density fluctuations, we distribute a pre- 
scribed number of particles A^o over two boxes, and de- 
vise an operation that exchanges particles between these 
boxes to maintain chemical equilibrium. By adopting 
the symmetrical restricted Gibbs (RG) ensemble 10], in 
which the boxes have equal constant volumes V = L'^, 
geometric cluster moves can be used for this exchange 
operation. The prescription for a cluster move closely 
follows the original algorithm [3, with the crucial dif- 
ference that a geometric operation not only alters the 
position of a particle, but also transfers it from one box 
to the other. Specifically, a pivot is placed at a ran- 
dom position within the first box, as well as at the corre- 
sponding position within the second box. A particle i is 
picked at random from one of the boxes (denoted 1) and 
point-reflected with respect to the pivot from its original 
position ri to r^. However, instead of placing the parti- 
cle at , we place it at the corresponding point in the 
other box (denoted 2), subject to periodic boundary con- 
ditions. Thereafter, any particle j interacting with par- 
ticle i around its original position in box 1 or its new po- 
sition in box 2 will also be considered for point reflection 
around the pivot and subsequent transfer to the opposite 
box, with a probability pij = max[l — exp(— /3Aij), 0], 
where Aij = —V{\ri — rj\) if particles i and j originally 
reside in the same box and Ay = V{\f[ — rj\) if i and j 
originally reside in different boxes. V{r) denotes a gen- 
eral pair potential and f3 the inverse temperature l/fceT. 
Note that pij solely depends on the pair interaction be- 
tween particles i and j, rather than on the total energy 
change resulting from the displacement of particle j . The 
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cluster construction proceeds iteratively for all particles 
interacting with each particle j. Upon completion of a 
cluster move, a new pivot is selected. The proof of de- 
tailed balance is analogous to that for the generalized 
GCA 0,0. 

The exchange of particles between boxes leads to 
(complementary) density fluctuations around the aver- 
age number density po = No/{2V) in each box. The 
fluctuation spectrum of the number density in box 1, 
pi — Ni/V, can be related to the grand- canonical prob- 
ability distributions P of the number densities of both 
boxes via 

P''''{pi\po,V,T) ^ P{p,\^l,V,T)P{2po - pi\p,V,T) , 

(1) 

where we note that P^'-^{pi) is independent of the choice 
of chemical potential p on the right-hand side 

Since the right-hand side of Eq. (Q) is the product of 
a function and a shifted (by 2po) and reflected form of 
this function, P^'^(pi) is symmetric (even) with respect 
to its mean pi = po. To facilitate comparison of the 
forms of P^-^{pi) for various choices of po, it is expedient 
to consider distributions of zero mean, to which end we 
define x = pi — po and write 

P'^'^ix\po,T) ^ Pix + po\T)P{-x + po\T) , (2) 

where we have suppressed reference to the constant vol- 
ume V and the (arbitrary) chemical potential p. 

The parameter space of the RG ensemble is spanned 
by Po and T. In the vicinity of the critical point (pg,T'^) 
terminating a line of liquid-vapor coexistence of a pure 
fluid or fluid mixture, P^<^ exhibits universal scaling be- 
havior. Introducing reduced variables go = {po ^ Po)/Po 
and t = (T — T'^)/T'^, we make the following ansatz for 
the finite-size scaling (ESS) properties of P^'^(x), 

P^^{x\go,t) ^ aoLl^'''P^^{aoL'^'''x,axQoL^''' ,a2tL^/'') . 

(3) 

Here P^'^ is a universal function which is symmetric in 
X for all values of t and go] P and v are critical expo- 
nents, and flo, ai, 02 are nonuniversal metric factors. The 
arguments in t and go control deviations from critical- 
ity. The temperature field has the form familiar from 
the ESS properties of magnets |0| or fluids 0, while 
that in po is particular to the RG ensemble. As one can 
verify 141 from an expansion of Eq. ^ with respect to 
Po, together with the known [l^ symmetry properties of 
the derivatives of P{x), variations in the form of P^'^{x) 
are — to leading order — controlled by the value of p^\ all 
terms having odd powers of po are antisymmetric in x 
and hence absent on symmetry grounds. 

To characterize the form of P^^{x) as a function of go 
and t, it is useful to consider the behavior of the dimen- 
sionless fourth-order cumulant ratio Q = (a;^)^/(a;^) [l^ . 
whose scaling properties follow from Eq. JSJ as 



with Q a universal function and Q{0,0) = Q* . The value 
of Q* = 0.711901 is known a priori by virtue of the re- 
suh that P^^{aoLf^/''x,0,0) oc [P*(L^/'^m)]^ HI, with 
P*{L^/^rn) the universal critical Ising magnetization dis- 
tribution 15]. Measurements of QLigo,t) for a range of 
global densities po provide a useful route to locating crit- 
icality. Specifically, consider the locus of points in go-t 
space for which Qhigo, t) = Q* , which we term the "iso- 
Q* curve." Expanding Eq. Q with respect to t < and 
go, and recalling that only terms involving even powers 
of go are nonzero, one has 

QUgo, t) = Q* [1 + qiglL^^/'' + gaii'/"^ + 0{gt, t^) 

from which it follows that, sufficiently close to the critical 
point, the iso-Q* curve is a parabola in go-t space. 



gl^-{qML('-'(^y'^t 



(6) 



The maximum of this parabola (at go = t ^ 0) coin- 
cides with the critical point and hence, by fitting to a 
few estimated points on the iso-Q* curve, one can read- 
ily determine the critical parameters and T'^. 

Turning now to the task of obtaining subcritical coex- 
istence properties within the RG framework, we consider 
the peak positions of P^'~^{pi\po, t) on some (subcritical) 
isotherm. One finds that when po equals the coexistence 
diameter pd = {pg -\- pi)/2, with pg and pi the gas and liq- 
uid densities, the peak positions of P^'^[pi) coincide with 
the coexisting densities, which can thus be simply read off 
from a measured histogram of its form An effective 
method for locating the diameter pd exploits the fact that 
the even moments of P^'^{x\po,t) are maximized when 
Po = Pd- From the absence of odd powers of po in the 
expansion of Eq. Q , it can then be shown that the 
variance a'^{po\t) of P^'^{x\po,t) varies to leading order 
quadratically in {po — Pd), i-e., 



'^^(PolO = <7^{Pd,t) - b{po - pdf 



(7) 



Ago,t)^Q{qigoLP^'',q2tL'^n, 



(4) 



with b a positive constant. By fitting to estimates of 
a'^{po\t), this result facilitates a determination of pd and 
thence the coexisting densities. 

To test the scaling theory, we first perform simulations 
using the GCA in the restricted Gibbs ensemble for a 
pure Lennard- Jones (LJ) fluid. We employ a potential 
cutoff 2.5(7 and reduced system sizes L* = L/a = 10, 
20, 30. At flxed global density po, histogram reweight- 
ing can be used to determine a point on the iso-Q* 
curve, i.e., the temperature at which Ql{po,T) takes 
the value Q*. Repeating for a range of po values al- 
lows the entire iso-Q* line to be mapped. As shown in 
Fig. n] the iso-Q* curves for the various L* are indeed 
parabolas [cf. Eq. lIHll] that coincide at their maximum. 
A careful analysis |l4| of the position of this maximum, 
f" = k^T^/e = 1.1878 (2), p" = pgcr^ = 0.3204 (5), re- 
veals excellent agreement with an existing GC estimate, 
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FIG. 1: (color online) Measured iso-Q* points (symbols) for 
a pure LJ fluid, together with parabolic fits (curves). The 
maxima of the curves (shown for L/cr = 10, 20, 30) in the po- 
T plane correspond to the critical point. The data collapse in 
the inset confirms the scaling prediction [Eq. 0]. 



= 1.1876 (3), = 0.3197 (4) PJ]. The inset confirms 
the finite-size scaling predicted in Eq. © with remark- 
able accuracy. Furthermore, the critical density distri- 
bution P^^{x) is indeed in quantitative agreement (not 
shown) with the square of the critical Ising magnetization 
distribution, as predicted pTl |. 

Having established the validity of our methodology, we 
now exploit the strengths of the GCA to address a typ- 
ical problem that is intractable for conventional simula- 
tions in the GC, NpT, or Gibbs ensemble. We consider 
a binary, strongly size-asymmetric mixture of LJ parti- 
cles of size a and small particles ( "additives" ) of diameter 
CTs = a/10. Depending on their interaction with the large 
particles, the presence of additives will affect the phase 
behavior and shift the location of the liquid- vapor critical 
point compared to the pure fluid. The additives mutually 
interact via a weakened LJ potential. 



(8) 

whereas a large and a small particle interact as hard 
spheres at a separation ct/j = {a + o's)/2. 

As before, we perform simulations for a range of po, 
at fixed additive volume fraction ips = ^cr^Ps = 0.005, 
corresponding to Ns = 10000. Because the small par- 
ticles are so numerous, and disperse relatively homoge- 
neously, the insertion probability of a fluid (large) par- 
ticle in a standard GC approach would be prohibitively 
small. By contrast, the present scheme renders it fea- 
sible to equilibrate the system and sample the density 
fluctuations. Figure |2] (diamonds) shows that the iso-Q* 
curve (plotted as a function of the total reduced density 
Po = pia^ + PsCTg) has a parabolic shape, as for the pure 
fluid. However, despite the relatively small volume frac- 
tion of additives, the maximum of this curve, i.e., the 
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FIG. 2: (color online) Iso-Q* curves for size-asymmetric bi- 
nary mixtures consisting of a LJ fluid (particle size a) and 
small additives (particle diameter cr/lO, volume fraction (jts)- 
All curves pertain to a linear system size L = 8.06cr. For addi- 
tives that interact with the fluid as hard spheres (diamonds, 
upper curve), the critical temperature and density increase 
compared to the pure fluid (open circles) , whereas the critical 
temperature decreases strongly if the additives have a weak 
attraction with the fluid (squares, filled circles, triangles). 



liquid- vapor critical point of the mixture, is markedly 
shifted. The increase in the critical temperature reflects 
an enhanced attraction between the large particles which 
stems from the entropic depletion interactions induced by 
the additives (isf. 

To highlight the subtle role of the nature of the inter- 
actions of the additives with the larger species, as well 
as their concentration, we study several systems in which 
there is a weak attraction between large and small par- 
ticles. The interaction is again of the LJ form, Eq. ||SJ), 
in which as is replaced with ais- Already at (jig = 10^^ 
the critical temperature is now noticeably depressed, as is 
evident from the shift of the iso-Q* maximum in Fig. |21 
and at 0s = 10~^, T'^ has decreased by almost 20%! 
We explain this surprising effect by the formation of a 
shell of small particles around the large particles, akin to 
nanoparticle halo formation |3i S j which weakens the 
effective attraction between the large particles. 

Our approach not only yields accurate estimates of 
critical points, but also entire coexistence curves. As 
described above, for each subcritical temperature, the 
variance of P^^{x) has a maximum at the coexistence 
curve diameter pd [Eq. Q], as is confirmed in Fig. |31 
The (total) densities of the coexisting liquid and vapor 
phases are determined from the peak positions of P^'^ 
[see Fig. EJa) for an example] , resulting in the phase di- 
agram in Fig.^Jb). We emphasize that obtaining such a 
phase diagram in a reasonable timescale would not be fea- 
sible using even the most efficient traditional approach to 
fluid phase equilibria, namely GC simulation 17]. Our 
tests show that the GC relaxation time is too large to 
be reliably estimated. Nevertheless, a lower bound on 
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FIG. 3: (color online) Variance of P^^{x) for a LJ mixture 
with size asymmetry a /as — 10, as a function of total number 
density po- Each (isothermal) curve is obtained from simula- 
tions at selected densities, but with constant additive volume 
fraction (j>3 = 0.005, and is parabolic (see fits), confirming 
Eq. 0. The maxima locate the coexistence curve diameter. 
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FIG. 4: (a) Example of a at po = Pd, for the fluid mixture 
described in Fig. |3 (b) Corresponding phase diagram; the 
open circle indicates the critical point of a pure LJ fluid. The 
full line is a fit of the form p^ — p" — ut + vt^ . 



the GC relaxation time, relative to that of the pure LJ 
fluid, can be estimated via a comparison of the large- 
particle transfer (insertion/deletion) acceptance proba- 
bility Pace- For liquid-hke densities of the large particles 
[pi « 0.6), we find that for (/)^ = 0.005, Pacc lO"'^; 
while for 0^ — 0.01, this falls to Pacc ^ 10"^. These val- 
ues are to be compared with pacc ~ 10~^ for the pure LJ 
fluid. One can therefore expect the GC relaxation time 
of the mixtures we have studied to be several orders of 
magnitude greater than for the pure LJ fluid. Since the 
algorithm presented here operates along fundamentally 
distinct lines — proposing large-scale collective updates of 
clusters of small and large particles, and accepting them 
with unit probability — it is not hampered by this prob- 
lem. Consequently it allows the efficient calculation of 
phase diagrams, even under conditions for which the GC 
approach fails. 



Summarizing, we have extended the rejection-free 
GCA to the study of phase transitions, by embedding it 
within a restricted Gibbs ensemble. The accurate loca- 
tion of critical points and coexistence curves within this 
ensemble requires a suitable FSS theory, which has been 
presented as well. By means of illustration, we have ap- 
plied our method to a strongly size-asymmetric LJ mix- 
ture, which cannot be studied with existing direct meth- 
ods, e.g., the GC ensemble. We find that the liquid- vapor 
phase behavior is highly sensitive to the concentration of 
small particles and the nature of their interaction with 
the large ones. Thus our method should prove useful in 
predicting the alterations to phase behavior which occur 
when small particles of various types are added to a fluid. 
Furthermore, by employing the method with a variant of 
the GCA suitable for electrostatic interactions T^, it be- 
comes possible to study the effects of adding salt on the 
phase behavior of charged colloids. Finally, we note that 
while for concreteness we have developed the formalism 
for the case of phase transitions whose order parameter 
is the density, the structure of our theory holds also for 
consolute points, or indeed situations where the order 
parameter is a linear combination of density and concen- 
tration. 

This material is based upon work supported by the 
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